
function [fPlus, fCross, dFPlusOverdTheta, dFPlusOverdPhi, dFCrossOverdTheta, dFCrossOverdPhi] = CalculateFpFcAndDerivatives(phi, theta, psi, Det)

%% This function calculates the antenna pattern functions Fp and Fc for any detector at a given sky location (phi, theta) and 
%% polarization angle (psi). The Det input is obtained by a call to LoadDetectorData in matapps which populates the structure.
%% A typical call to this function would be :
%% Det = LoadDetectorData('INDIGO')
%% [fp, fc, dfpdTheta, dfpdPhi, dfcdTheta, dfcdPhi] = CalculateFpFcAndDerivatives(phi, theta, psi, Det);
%% Note that this function is not vectorized - and can handle only one direction (phi, theta) at a time. 


 d11 = Det.d(1,1); 
 d12 = Det.d(1,2);
 d13 = Det.d(1,3); 
 d21 = Det.d(2,1); 
 d22 = Det.d(2,2); 
 d23 = Det.d(2,3); 
 d31 = Det.d(3,1); 
 d32 = Det.d(3,2); 
 d33 = Det.d(3,3); 

 fPlus = ...
   d11*(Power(cos(psi)*sin(phi) - cos(phi)*cos(theta)*sin(psi),2) - Power(cos(phi)*cos(psi)*cos(theta) + sin(phi)*sin(psi),2)) + ... 
   d12*(-(cos(2*psi)*(3 + cos(2*theta))*sin(2*phi))/4. + cos(2*phi)*cos(theta)*sin(2*psi)) + ... 
   d21*(-(cos(2*psi)*(3 + cos(2*theta))*sin(2*phi))/4. + cos(2*phi)*cos(theta)*sin(2*psi)) + ...
   d22*(Power(cos(phi),2)*cos(2*psi) + cos(theta)* (-(cos(2*psi)*cos(theta)*Power(sin(phi),2)) + sin(2*phi)*sin(2*psi))) + ...
   d23*(Power(cos(psi),2)*cos(theta)*sin(phi) - 2*cos(phi)*cos(psi)*sin(psi) - cos(theta)*sin(phi)*Power(sin(psi),2))*sin(theta) + ...
   d32*(Power(cos(psi),2)*cos(theta)*sin(phi) - 2*cos(phi)*cos(psi)*sin(psi) - cos(theta)*sin(phi)*Power(sin(psi),2))*sin(theta) + ... 
   d13*(cos(phi)*cos(2*psi)*cos(theta) + sin(phi)*sin(2*psi))*sin(theta) + ...
   d31*(cos(phi)*cos(2*psi)*cos(theta) + sin(phi)*sin(2*psi))*sin(theta) - ...
   d33*cos(2*psi)*Power(sin(theta),2);
     

 fCross = ...
  -2*d11*(cos(psi)*sin(phi) - cos(phi)*cos(theta)*sin(psi))*(cos(phi)*cos(psi)*cos(theta) + sin(phi)*sin(psi)) + ...
   2*d22*(cos(psi)*cos(theta)*sin(phi) - cos(phi)*sin(psi))*(cos(phi)*cos(psi) + cos(theta)*sin(phi)*sin(psi)) + ...
     d12*(Power(cos(phi),2)*cos(2*psi)*cos(theta) - cos(2*psi)*cos(theta)*Power(sin(phi),2) + ((3 + cos(2*theta))*sin(2*phi)*sin(2*psi))/4.) + ...
     d21*(Power(cos(phi),2)*cos(2*psi)*cos(theta) - cos(2*psi)*cos(theta)*Power(sin(phi),2) + ((3 + cos(2*theta))*sin(2*phi)*sin(2*psi))/4.) + ... 
     d13*(Power(cos(psi),2)*sin(phi) - 2*cos(phi)*cos(psi)*cos(theta)*sin(psi) - sin(phi)*Power(sin(psi),2))* sin(theta) + ...
     d31*(Power(cos(psi),2)*sin(phi) - 2*cos(phi)*cos(psi)*cos(theta)*sin(psi) - sin(phi)*Power(sin(psi),2))*sin(theta) - ...
     d23*(cos(phi)*cos(2*psi) + cos(theta)*sin(phi)*sin(2*psi))*sin(theta) - ...
     d32*(cos(phi)*cos(2*psi) + cos(theta)*sin(phi)*sin(2*psi))*sin(theta) + ...
     d33*sin(2*psi)*Power(sin(theta),2);

 dFPlusOverdTheta = ...
   (d23 + d32)*Power(cos(psi),2)*cos(2*theta)*sin(phi) - d23*Power(cos(theta),2)*sin(phi)*Power(sin(psi),2) - ... 
    d32*Power(cos(theta),2)*sin(phi)*Power(sin(psi),2) + d13*cos(theta)*sin(phi)*sin(2*psi) + ... 
    d31*cos(theta)*sin(phi)*sin(2*psi) + cos(phi)* ((d13 + d31)*cos(2*psi)*cos(2*theta) - (d23 + d32)*cos(theta)*sin(2*psi)) - ... 
    d12*cos(2*phi)*sin(2*psi)*sin(theta) - d21*cos(2*phi)*sin(2*psi)*sin(theta) + ... 
    d11*sin(2*phi)*sin(2*psi)*sin(theta) - d22*sin(2*phi)*sin(2*psi)*sin(theta) + ...
    d23*sin(phi)*Power(sin(psi),2)*Power(sin(theta),2) + d32*sin(phi)*Power(sin(psi),2)*Power(sin(theta),2) - ...
    d33*cos(2*psi)*sin(2*theta) + d11*Power(cos(phi),2)*cos(2*psi)*sin(2*theta) + ...
    d22*cos(2*psi)*Power(sin(phi),2)*sin(2*theta) + (d12*cos(2*psi)*sin(2*phi)*sin(2*theta))/2. + ...
    (d21*cos(2*psi)*sin(2*phi)*sin(2*theta))/2. ;

 dFPlusOverdPhi = ...
   (-3*d22*cos(2*psi)*sin(2*phi))/2. - (d22*cos(2*psi)*cos(2*theta)*sin(2*phi))/2. - ...
   4*d11*Power(cos(phi),2)*cos(psi)*cos(theta)*sin(psi) + 4*d11*cos(psi)*cos(theta)*Power(sin(phi),2)*sin(psi) - ... 
   2*d12*cos(theta)*sin(2*phi)*sin(2*psi) - 2*d21*cos(theta)*sin(2*phi)*sin(2*psi) - ... 
   (cos(2*phi)*((d12 + d21)*cos(2*psi)*(3 + cos(2*theta)) - 4*d22*cos(theta)*sin(2*psi)))/2. - ...
   d13*cos(2*psi)*cos(theta)*sin(phi)*sin(theta) - d31*cos(2*psi)*cos(theta)*sin(phi)*sin(theta) + ... 
   d23*sin(phi)*sin(2*psi)*sin(theta) + d32*sin(phi)*sin(2*psi)*sin(theta) + ... 
   cos(phi)*((d13 + d31)*sin(2*psi)*sin(theta) + cos(2*psi)*(d11*(3 + cos(2*theta))*sin(phi) + (d23 + d32)*cos(theta)*sin(theta))) ;

 dFCrossOverdTheta = ...
  (d13 + d31)*Power(cos(psi),2)*cos(theta)*sin(phi) - d13*cos(theta)*sin(phi)*Power(sin(psi),2) - ...
   d31*cos(theta)*sin(phi)*Power(sin(psi),2) - d23*Power(cos(theta),2)*sin(phi)*sin(2*psi) - ... 
   d32*Power(cos(theta),2)*sin(phi)*sin(2*psi) + d12*cos(2*psi)*Power(sin(phi),2)*sin(theta) + ...
   d21*cos(2*psi)*Power(sin(phi),2)*sin(theta) - 4*d22*cos(psi)*cos(theta)*Power(sin(phi),2)*sin(psi)*sin(theta) - ...
   Power(cos(phi),2)*((d12 + d21)*cos(2*psi) + 4*d11*cos(psi)*cos(theta)*sin(psi))*sin(theta) + ...
   d23*sin(phi)*sin(2*psi)*Power(sin(theta),2) + d32*sin(phi)*sin(2*psi)*Power(sin(theta),2) + ...
   d33*sin(2*psi)*sin(2*theta) - cos(phi)*((d23 + d32)*cos(2*psi)*cos(theta) + ...
   (d13 + d31)*Power(cos(theta),2)*sin(2*psi) - 2*d11*Power(cos(psi),2)*sin(phi)*sin(theta) + ...
   2*d22*Power(cos(psi),2)*sin(phi)*sin(theta) + 2*d11*sin(phi)*Power(sin(psi),2)*sin(theta) - ...
   2*d22*sin(phi)*Power(sin(psi),2)*sin(theta) - d13*sin(2*psi)*Power(sin(theta),2) - ...
   d31*sin(2*psi)*Power(sin(theta),2) + d12*sin(phi)*sin(2*psi)*sin(2*theta) + ...
   d21*sin(phi)*sin(2*psi)*sin(2*theta));

dFCrossOverdPhi = ...
  -2*(d11 - d22)*Power(cos(phi),2)*cos(2*psi)*cos(theta) +  ...
   cos(phi)*((d13 + d31)*cos(2*psi)*sin(theta) - 2*cos(psi)*sin(psi)*((d11 - d22)*(3 + cos(2*theta))*sin(phi) + (d23 + d32)*cos(theta)*sin(theta))) + ... 
   (4*(d11 - d22)*Power(cos(psi),2)*cos(theta)*Power(sin(phi),2) - ...
      4*d11*cos(theta)*Power(sin(phi),2)*Power(sin(psi),2) + ...
      4*d22*cos(theta)*Power(sin(phi),2)*Power(sin(psi),2) + 3*d12*cos(2*phi)*sin(2*psi) + ...
      3*d21*cos(2*phi)*sin(2*psi) + d12*cos(2*phi)*cos(2*theta)*sin(2*psi) + ...
      d21*cos(2*phi)*cos(2*theta)*sin(2*psi) + ...
      2*cos(2*psi)*sin(phi)*(-4*(d12 + d21)*cos(phi)*cos(theta) + (d23 + d32)*sin(theta)) + ...
      d13*sin(phi)*sin(2*psi)*sin(2*theta) + d31*sin(phi)*sin(2*psi)*sin(2*theta))/2.;


